library(plotly)
## Warning: package 'plotly' was built under R version 3.3.3
## Loading required package: ggplot2
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
library(microbenchmark)
## Warning: package 'microbenchmark' was built under R version 3.3.3
grid <- read.csv("~/GitHub/SDS385-course-work/Excercise 8 Spatial smoothing/Laplacian smoothing/fmri_z.csv")
grid=as.matrix(grid)
source('~/GitHub/SDS385-course-work/Excercise 8 Spatial smoothing/Laplacian smoothing/functions.R')
p1=plot_ly(z=grid,type='heatmap')
n=nrow(grid)
y=as(c(grid),"CsparseMatrix")
lambda=0.1
C=makeC_sparse(n,n,lambda)
x = solve(C,y,sparse=TRUE)
smoothed_grid = matrix(data=x,nrow=n)
p2=plot_ly(z=smoothed_grid,type='heatmap')
lambda=1
C=makeC_sparse(n,n,lambda)
x = solve(C,y,sparse=TRUE)
smoothed_grid = matrix(data=x,nrow=n)
p3=plot_ly(z=smoothed_grid,type='heatmap')
lambda=10
C=makeC_sparse(n,n,lambda)
x = solve(C,y,sparse=TRUE)
smoothed_grid = matrix(data=x,nrow=n)
p4=plot_ly(z=smoothed_grid,type='heatmap')
p=subplot(p1,p2,p3,p4,nrows=2)
p
lambda=1
C=makeC_sparse(n,n,lambda)
result=gauss_seidel(C,y,50,targetfunction = TRUE)
targetfunction=result[[2]]
plot(targetfunction,type='l',col='blue')
legend("topright",col=c('blue','red','black'),legend=c("Gauss Seidel","Jacobi","Direct solver"),lty=1)
result=jacobi(C,y,50,targetfunction = TRUE)
targetfunction=result[[2]]
lines(targetfunction,col='red')
x = solve(C,y,sparse=TRUE)
direct_solver=target_function(C,y,x)[1,1]
lines(rep(direct_solver,50),col='black')

x=gauss_seidel(C,y,8)[[1]]
smoothed = matrix(data=x,nrow=n)
plot_ly(z=smoothed,type='heatmap')
x=jacobi(C,y,15)[[1]]
smoothed = matrix(data=x,nrow=n)
plot_ly(z=smoothed,type='heatmap')
microbenchmark(gauss_seidel(C,y,8),times=1)$time
## [1] 28099008
microbenchmark(jacobi(C,y,15),times=1)$time
## [1] 31740398
microbenchmark(solve(C,y,sparse=TRUE),times=1)$time
## [1] 4680716